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Section  1 


INTRODUCTION 


The  application  of  linear  elastic  fracture  mechanics  analysis 
to  structural  details  in  new  aircraft  designs  has  received  growing 
emphasis  from  the  Air  Force  and  from  aircraft  manufacturers  during 
the  past  few  years.  Not  only  have  fracture  mechanics  data  become 
more  readily  available  in  recent  years  [1,2] , but  also,  there  has 
been  a trend  toward  treatment  of  the  problem  of  fracture  in  its 
own  right,  distinct  from  fatigue.  Within  the  past  year,  this 
trend  has  culminated  in  the  establishment  of  structural  design 
criteria  by  the  Air  Force  which  require  an  aircraft  designer  to 
take  specific  design-analysis  actions  to  protect  structures  against 
fracture  [3] . Required  design  calculations  now  include,  generally, 
comparison  of  load-induced  stress  intensity  factors  to  material 
fracture  toughness  and  assessment  of  crack  growth  rates,  based  on 
assumptions  concerning  the  size  and  location  of  possible  cracks  in 
the  structure.  Both  types  of  calculation  require  prior  estimation 
of  the  load-induced  stress  intensity  factor.  Hence,  there  has  also 
been  considerable  emphasis  on  adding  to  the  body  of  available  fracture 
mechanics  solutions. 

Because  so  few  geometrical  configurations  are  amenable  to  a 
purely  analytical  solution  of  the  equations  of  elastic  fracture 
mechanics,  the  finding  of  new  solutions  depends  upon  development 
of  numerical  analysis  techniques.  Extensive  contributions  have 
been  made  by  Bowie  and  his  colleagues  [4,5,6]  using  the  complex 
variable  formulation  of  elasticity  in  combination  with  conformal 
mapping,  analytic  continuation  and  boundary  collocation  methods. 

Tada,  et  al.  [7]  have  recently  collected  and  classified  a compre- 
hensive body  of  solutions  based  on  the  semi-analytical  methods 
(complex  variables,  boundary  collocation,  Fourier  transforms,  etc.) 
among  which  appear  many  new  solutions  by  Tada.  However,  the  semi- 
analytical  methods  have  not  as  yet  proved  capable  of  application 
to  the  irregular  geometrical  configurations  which  are  found  so  often 
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in  real  airframe  structural  details.  Stress  analysis  of  irregular 
structure  has  been  the  province  of  the  f inite-element  methods  for 
the  past  decade.  Within  the  last  six  years,  numerous  contributors 
have  extended  the  finite-element  technique  to  fracture  mechanics 
analysis . 

Work  on  finite-element  fracture  mechanics  analysis  at  MIT  has 
followed  the  path  of  assumed-stress  hybrid  elements,  first  proposed 
by  Pian  [8]  for  ordinary  continuum  elements.  The  hybrid  method  was 
subsequently  extended  by  Tong,  Pian,  Luk,  and  Lasry  [9,10,11,12]  to 
formulation  of  rectangular  elements  which  incorporate  an  elastic 
crack-tip  singularity,  but  which  also  have  an  assumed  linear  or 
quadratic  variation  of  displacements  along  the  element  edges  and 
are  thus  compatible  with  ordinary  elements  (Figure  1) . Numerical 
experiments  have  demonstrated  that  the  hybrid  crack-containing  ele- 
ments are  capable  of  producing  estimates  of  and  K with  less 
than  1 percent  error,  using  20  to  50  total  degrees  of  freedom  in 
the  analysis,  for  simple  geometrical  configurations.  Hence,  it 
becomes  possible  to  create  economically  practical  analysis  procedures 
for  structural  details  by  refining  the  mesh  or  ordinary  elements  to 
pick  up  the  stress  gradients  caused  by  nonuniform  loading  and  com- 
plicated boundary  geometry,  leaving  to  the  special  hybrid  element 
the  task  of  picking  up  the  local  gradients  caused  by  the  crack-tip 
singularity. 

This  report  summarizes  recent  developments  at  the  MIT  Aero- 
elastic  and  Structures  Research  Laboratory  (ASRL)  in  which  the  crack- 
containing  hybrid  element  has  been  applied  for  the  first  time  to  some 
typical  structural  details,  found  in  current  production  aircraft, 
with  geometries  too  complicated  for  economical  solution  by  other 
techniques.  The  "PCRK59"  crack  element  used  in  these  analyses  is 
a generalized  version  of  the  original  Lasry  element  [12]  which  was 
formulated  and  programmed  by  Tong  and  subsequently  modified  for 
greater  utility  by  the  ASRL  computing  staff. 


2- 


Section  2 


BASIC  ELEMENTS  AND  METHODS 


2 . 1 Element  QUAD4 

The  ASRL  QUAD4  four-node  quadrilateral  element  (Figure  2)  is 
used  as  the  basic  building  block  in  the  analysis  procedure.  QUAD4 
is  the  well-known  bilinear  isoparametric  assumed-displacement  ele- 
ment which  has  been  used  for  continuum  stress  analysis  for  many 
years  [13] . The  ASRL  version  has  been  programmed  as  an  independent 
subroutine  which  includes  the  options  of  individual  rotation  trans- 
formations at  each  node  and  calculation  of  a "B"  matrix  for  stress 
analysis. 

The  nodal  coordinates  X^ , , X2,...,Y^,  element  thickness  T 

and  the  elastic  constants  matrix: 
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Subroutine  QUAD4  allows  for  general  plane  orthotropic  behavior,  a 
capability  which  is  included  in  the  attachment  lug  program.  Plane 
stress  is  assumed  in  the  analysis.  The  lug  program  does  not  use 
the  rotation  transformation  option. 

The  element  stiffness  matrix  k is  calculated  by  numerical 
area  integration,  using  3x3-point  Gaussian  quadrature  [14]  for 
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ELEMENT 


l<=  t //  2T£i  Jxjv  ,3) 

where  D contains  the  interior  strain-nodal  displacement  relations: 

fox  Err  %z  •/g'J  (4) 

The  stiffnesses  are  returned  in  Lower  Triangle  Vector  (LTV)  form: 

l<  = L *<n  k3l  k22  k31  kJZ  ■ ...  k8ff  j (5) 

For  the  purpose  of  stress  analysis,  Eqs.  1 and  4 may  be  combined 
to  give: 

r(x,Y).  cj(x,y)f=  B(x,y>j.  (6) 


QUAD4  also  returns  the  matrix  B(X  ,Y  ),  formed  at  the  fifth 

- c c 

Gaussian  station,  for  later  calculation  of  stresses  at  the  element 
"centroid",  defined  in  terms  of  the  nodal  coordinates: 


Xc 


4 

= x r 

+ i=i 


Yc  * |:L  Y; 
C '1=1 


(7) 


The  behavior  of  QUAD4  has  been  studied  extensively  on  other 
projects  and  is  well  understood.  Uniform  or  nearly  uniform  stress 
fields  can  be  picked  up  to  within  the  roundoff  accuracy  of  the  digi- 
tal computer  being  used  for  the  analysis.  The  inability  of  the 
bilinear  assumed  displacement  fields  to  follow  the  quadratic  deflec- 
tion of  the  neutral  axis  of  a cantilever  beam  loaded  by  an  end  moment 
has  been  well  documented  elsewhere  [15]  and  constitutes  a limitation 
on  the  QUAD 4 . In  practical  terms,  this  requires  that  the  element 
aspect  ratio  (Figure  2)  be  held  close  to  unity  for  models  of  struc- 
tures which  are  expected  to  have  quadratic  or  higher-order  displace- 
ment behavior.  In  some  cases,  even  an  aspect  ratio  of  unity  is  not 
sufficient  to  insure  convergence  of  the  solution.  For  example,  the 
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QUAD 4 element  was  used  recently  to  model  a thick-walled  cylinder 

subjected  to  centrifugal  loading  from  its  own  mass,  due  to  rotation 

[16] . The  analytical  solution  of  this  axisymmetric  problem  includes 
3 

an  r term  in  the  radial  displacement  field  which  the  bilinear  ele- 
ment is  unable  to  pick  up;  errors  of  25%  were  found  for  a cylinder 
with  a 2:1  ratio  of  outside-to-inside  radius,  using  four  unit-aspect- 
ratio  QUAD 4 elements  through  the  wall  thickness. 

The  misbehavior  of  the  bilinear  element  in  the  presence  of 
higher-order  gradients  requires  the  use  of  many  elements  to  model 
complicated  geometries.  Also,  "calibration"  of  the  finite-element 
model  is  a good  idea,  where  possible,  by  comparing  the  numerical 
results  with  independent  solutions.  Calibrations  for  this  project 
have  included  comparisons  with  the  classical  elasticity  solution 
for  stresses  and  displacements  near  a circular  hole  in  a semi- 
infinite strip  under  tension  [17]  and  with  finite-element  analyses 
using  higher-order  assumed-displacement  elements. 

2.2  Element  PCRK59 

Formulation  of  the  assumed-stress  hybrid  finite-element  method 
begins  with  the  Principle  of  Minimum  Complementary  Energy: 

Tfc  • £ [ f “TI  JS  - JifS  z M ] (8) 

" su  V 

where 

£ indicates  summation  over  the  element  set. 
n 

S = part  of  the  element  boundary  over  which  displace- 
ments are  prescribed. 

V = element  volume. 

u = vector  of  prescribed  displacements  on  S^. 
a = stress  vector. 

S = compliance  constants  = C ^ (e  = So) . 

T = vector  of  surface  tractions  = Na,  where  N is  a 
matrix  of  surface  normal  direction  cosines. 
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If  II  is  used  directly,  only  the  stress  field  is  assumed,  subject 
to  admissibility  criteria  requiring  that  the  assumed  stresses 
satisfy: 

(i)  Interior  equilibrium  = 0,  in  V, 

where  F^  are  prescribed  body  forces. 

A 

(ii)  Mechanical  boundary  conditions  Na  = T on  S , 
that  part  of  the  element  boundary  over  which 

/v 

the  surface  tractions  T are  prescribed. 

(iii)  Equilibrium  of  surface  tractions  Na  across  the 

interelement  boundaries  S,  which  are  distinct 

from  S and  S . 
u a 

Formal  application  of  the  variational  calculus  to  Eq.  8 leads  to 
two  sets  of  Euler  equations: 

(iv)  Interior  compatibility.  So  = e in  V,  where 

e..  = 1/2  (3u./3x.  + 3u./3x.). 
ij  i J 3 i 

(v)  Displacement  boundary  conditions,  u = u on  S^. 

If  assumed  stress  functions  are  substituted  and  Eq . 8 is  integrated 
before  Hc  is  varied,  there  results  a linear  equation  system  in  which 
the  generalized  coordinates  to  be  solved  for  are  forces.  This 
approach  leads  to  a Matrix  Force  Method  analysis  which  brings  with 
it  the  programming  problem  of  systematic  identification  and  elimina- 
tion of  redundant  quantities. 

The  assumed-stress  hybrid  approach  avoids  the  complications  of 
force  redundancy  by  modifying  II  so  that  the  primary  unknowns  in  a 
finite-element  application  become  displacements  once  again.  The 
Principle  of  Minimum  Complementary  Energy  is  modified  by  addition 
of  Lagrange  multiplier  terms  [8,18]  which  change  admissibility 
criteria.  Specifically,  conditions  (ii)  and  (iii)  above  are  relaxed 
and  confition  (v)  is  enforced.  Under  the  new  principle,  stress  func- 
tions satisfying  only  the  interior  equilibrium  conditions  (i)  may  be 
assumed,  and  displacement  functions  which  satisfy  interelement  com- 
patibility and  conditions  (v)  must  now  be  assumed  as  well.  The 
modified  energy  principle  which  replaces  Eq.  8 is: 
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(9) 


7T  = E [ fjT5^S  - J is’SZM  - £ »TJJs] 

i n igv  V -V 

where  "3V"  represents  the  entire  element  boundary,  S + + S . 

To  convert  11^  into  a finite-element  formulation,  the  stress  vector 
a is  assumed  within  each  element  and  the  displacement  vector  u is 
assumed  on  the  boundary,  3V,  of  each  element: 


a = P(x,y,z)tf  u = L (x,  y , z)  q (10) 

where  P and  L are  matrices  of  interpolation  functions.  Vector  3 
contains  generalized  stress  coordinates,  while  q is  a vector  of 
nodal  displacements.  Matrices  P and  L are  assumed  independently, 
with  L defined  only  along  the  element  boundary  3V.  Substitution 
of  Eq.  10  into  Eq.  9 then  leads  to: 


Vi*  $UT£%-i/!Tt!/l  -jre] 


(11) 


where 


S=  I Js  h=/ptsp^v 

dv  y 

Q=  J L TJS  = Consis+ewl  nodal  -force  vec+or 

~ 5 ~ 

-V 


(12) 


and  where  No  = NPB  has  been  substituted  for  the  surface  traction 
vector  T in  the  expression  for  G. 

Direct  assembly  and  solution  of  the  equation  system  represented 
by  ni  is  possible,  but  results  in  a mixed  matrix  method,  with  both 
force  and  displacement  unknowns.  A more  versatile  formulation  is 
obtained  by  recognizing  that,  since  the  stresses  are  assumed  inde- 
pendently within  each  element,  P3  for  one  element  is  not  coupled 
with  any  other  elements.  Therefore,  the  unknowns  3 may  be  formally 
eliminated  by  applying  the  variational  calculus  to  n^: 

m a TT 

61^(3  for  only  one  element  varied)  = 63  { 1/93)  = 

6BT(Gq  - H3)  = 0 
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which  leads  to 


-1 

/£=££$•  U3) 

-v/ 

* 

Substitution  of  Eq . 13  back  into  Eq.  11  then  yields  : 

X-  <n> 

* Y)  ~ Z,  ~ 

The  alternate  expression  for  I L given  by  Eq.  14  represents  a 

1 T - 1 

pure  Matrix  Displacement  Method.  The  quantity  G H G can  be  recog- 
nized as  an  element  stiffness  matrix.  The  E notation  may  now  be 

n 

identified  as  a conventional  Matrix  Displacement  Method  structure 

model  assembly  procedure.  Hence,  structure  models  can  be  created 

by  assembling  conventional  and  hybrid  elements,  provided  only  that 

compatibility  across  the  interelement  boundaries  is  maintained  by 

proper  choices  for  the  assemed  displacement  fields.  The  versatility 

of  the  hybrid  method  lies  in  its  ability  to  provide  special-purpose 

elements,  for  restricted  regions,  which  may  be  coupled  into  a model 

containing  conventional  elements  in  the  remaining  regions  which  are 

free  of  singularities  or  other  unusual  behavior. 

The  original  hybrid  crack  elements  [9,10]  were  derived  from 

-1/2 

ni  by  assuming  a stress  field  containing  r terms,  where  r 

measures  radial  distance  from  the  crack  tip  and,  by  assuming  dis- 
placement fields  which  vary  linearly  from  one  node  to  the  next, 
along  the  element  boundary.  However,  subsequent  analysis  of  error 

sources  [19]  has  indicated  that  the  area  integration  required  for 

-1/2 

computation  of  H (see  Eqs.  12)  gives  poor  results  for  the  r 
terms  since  some  of  the  Gaussian  stations  are  close  to  the  crack 
tip.  This  situation  may  be  remedied  by  increasing  the  number  of 
Gaussian  stations,  but  the  computation  of  k then  becomes  too  costly. 
A better  approach,  used  for  the  second-generation  crack  elements 
[11,12]  has  been  followed  in  the  present  work.  The  energy  principle 
11^  may  be  further  modified  by  introducing  two  displacement  fields: 
u assumed  on  3V  and  u assumed  in  V.  There  now  arises  another  com- 
patibility condition,  u = u on  8V,  which  is  relaxed  by  the  Lagrange 

* -1  . 

Note  that  H and  H are  symmetric  matrices. 
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multiplier  method.  At  the  same  time,  the  condition  that  u must 
satisfy  the  interior  equilibrium  equations,  as  well  as  the  strain- 
displacement  relations,  is  enforced.  As  a result,  the  area  inte- 
gration is  converted  to  a boundary  integral  and  is  modified  to 
the  form: 

7T  = z:f  J JtuJls  - i y i(ttu+iTt)js-  (15) 

The  same  boundary  displacement  field  u can  be  used  for  both  II  ^ and 
II ^ - However,  the  interior  assumed  fields  in  IT ^ must  be  a complete 
elasticity  solution:  stresses  a and  displacements  u which  satisfy 

all  of  the  equations  of  elasticity,  with  T = Na  a derived  quantity. 
The  distributions  for  a and  u are  obtained  from  a complex  variable 
solution  of  the  equations  of  elasticity  near  a crack  tip  or  equiva- 
lently, by  solving  the  biharmonic  equation  for  an  Airy  stress  func- 
tion. Computation  of  the  element  stiffness  matrix  is  the  same  as 
for  ni#  except  that  H is  now  computed  by  a boundary  integral: 

h-  1 HckDV  J5  <16) 

9V 

where  A is  a matrix  of  shape  functions  corresponding  to  the  interior 
assumed  displacement  field,  u = Aq. 

The  principle  II ^ also  possesses  the  advantage  of  convenience 
for  treatment  of  arbitrary  shapes,  since  only  boundary  integration 
is  required.  Figure  3 illustrates  the  PCRK59  element  based  on  n„. 
Input  information  required  by  this  element  is  similar  to  the  infor- 
mation required  by  QUAD4 : 

(i)  Geometry:  global  coordinates  of  the  crack  tip 

X, , Y. ; global  coordinates  of  each  node  X. , Yn ,X_ , . . . Y_ . 
t t l ± z y 

(ii)  Material  properties:  the  shear  modulus  G = E/2(l+v) 

and  a second  constant  rt  = (3-v)/(l+v)  for  plane  stress 

(3-4v  for  plane  strain) 

PCRK59  is  programmed  only  for  isotropic  material  and  does  not 
incorporate  the  rotation  transformations  available  in  QUAD4 . In 
fact,  rotation  transformations  cannot  be  applied  once  the  PCRK59 
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stiffness  matrix  has  been  formed.  This  limitation  is  caused  by 
the  appearance  of  the  crack  tip  coordinates  in  the  numerical  inte- 
gration scheme,  but  the  restriction  does  not  affect  many  practical 
fracture  mechanics  problems.  The  numerical  integration  is  by 
five-point  Gaussian  quadrature  [14]  between  each  pair  of  nodes, 
except  that  the  crack  surfaces  are  skipped.  Omission  of  the  crack 

/s 

surfaces  is  justified  because  they  constitute  S , over  which  T = 0, 
and  because  the  derived  tractions  T satisfy  this  stress-free  con- 
dition at  least  in  an  average  sense. 

The  PCRK59  element  has  two  other  important  features.  First, 
unit  thickness  is  assumed.  Second,  a symmetric  "half-element" 
option  is  available,  under  which  nodes  1,5  and  the  crack  tip  are 
assumed  to  lie  on  a line  parallel  to  the  global  X-axis,  while  the 
element  and  applied  loading  are  assumed  to  be  symmetric  about  this 
line.  Under  these  conditions,  a half-model  of  a structure  may  be 
analyzed  to  obtain  Mode  I stress  intensity  solutions  only?  e.g., 
for  the  coupon  in  uniform  tension  with  edge  cracks,  shown  in  Figure 
4.  The  "half-element"  consists  of  nodes  1,2,... 5,  only,  with  node 
1 requiring  a roller  restraint  to  maintain  the  assumed  symmetry. 
Another  input  parameter  determines  which  option  is  executed: 

KEY  = 1 for  "half-element" 

2 for  full  element 

The  "half-element"  option  is  used  mainly  for  illustrative  examples 
and  performance  testing.  The  full  element  option  has  been  used 
exclusively  in  the  present  work. 

Element  PCRK59  computes  and  returns  a stiffness  matrix  in  LTV 
form  (see  Eq.  5)  for  either  the  10  degree-of-f reedom  "half-element" 
or  the  18  degree-of-f reedom  full  element.  In  addition,  a special 
B matrix  for  calculation  of  stress  intensities  is  returned.  Eq.  13, 
used  in  the  derivation  of  the  stiffness  matrix,  can  also  be  used  to 
compute  the  generalized  stress  coordinates  £5  after  the  element  nodal 
displacements  q have  been  obtained.  For  the  PCRK59, 

Mi  A A •••/*?  A At  •••  (17) 
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where  k^k^r  are  the  Mode  I,  Mode  II  stress  intensity  factors, 
defined  by: 


£ (e)  + ^ZfB)  + (fermS  <n  • • > 


YTF 

Oyy  = etc. 


(18) 


The  functions  f^,f2,  are  from  the  classical  crack  tip  solution, 
and  the  other  generalized  coordinates  $2' ^3' * * * '^18  rePresen,,: 
far-field  behavior.  Thus,  B is  formed  by  extracting  the  first  and 
tenth  rows  of  H ^G,  so  that  the  stress  intensities  may  be  calculated 
from: 


l yki  ^ J-f  CluiB)  ^ ^ 


(19) 


for  the  full  element.  Only  the  first  row  of  H "*"G  is  extracted  if 
the  "half-element"  option  is  in  effect: 


- 5(lXio)  ^ t-2 


(20) 


NASA/ASTM  standard  stress  intensity  factors  may  be  computed  after 
Eq.  19  or  Eq.  20  by: 

Kj  = -fa?  Ki  - “k z (2i) 

If  a structure  model  with  thickness  T ^ 1 is  to  be  analyzed,  this 
may  be  done  simply  by  scaling  the  PCRK59  stiffness  to: 

k'r  T k (22) 

Performance  of  the  PCRK59  element  has  been  tested  extensively 
by  comparison  with  classical  and  boundary  collocation  solutions 
[19].  Solutions  for  accurate  to  better  than  1 percent  have  been 
obtained  with  a rectangular  crack  element  surrounded  by  only  a few 
QUAD4  elements.  Other  tests  have  shown  that  solution  accuracy 
within  3 percent  is  maintained  when  the  crack  element  shape  is  dis- 
torted by  relocating  some  nodes  as  much  as  0.3  x (length  of  crack 
within  element)  away  from  the  positions  they  occupy  for  a rectangle. 
Also,  the  3 percent  accuracy  limit  can  be  maintained  with  the  crack 
tip  located  anywhere  from  20  to  70  percent  across  a line  between 
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nodes  5 and  1,  with  the  element  shape  kept  rectangular.  The  dis- 
tortions of  element  shape  and  crack  tip  location  required  for  the 
structure  models  analyzed  in  the  present  work  are  well  within  these 
limits . 

The  PCRK59  element  possesses  one  unavoidable  quirk  which  arises 
from  its  linearity.  If  the  element  is  placed  in  a region  with*  com- 
pressive stress  normal  to  the  crack,  a negative  value  of  K is 
obtained.  In  a real  structure,  the  crack  would  close  and  cease  to 
be  a problem  in  this  situation.  Therefore,  negative  K values 
should  be  interpreted  as  signaling  the  absence  of  Mode  I stress 
intensity.  On  the  other  hand,  the  solution  for  K will  be  posi- 
tive (negative)  according  to  whether  the  crack  is  being  subjected 
to  positive  (negative)  shear  stress,  as  defined  by  the  standard 
conventions  of  elasticity.  In  this  case,  the  correct  interpreta- 
tion is  to  take  the  absolute  value  of  K . 

In  summary,  the  PCRK59  element  permits  efficient  computation 
of  stress  intensity  factors  by  well  established  procedures  of  the 
Matrix  Displacement  Method.  The  unusual  features  of  the  element 
are  internal  to  its  subroutine.  The  element  subroutine  requires 
familiar  input  information  and  returns  k and  B matrices  like  a con- 
ventional element.  The  structure  model  is  assembled  and  a global 
displacement  solution  is  computed  by  standard  techniques.  Computa- 
tion of  either  the  centroid  stresses  in  the  conventional  elements 
or  the  stress  intensity  factors  in  the  crack  element  is  then  merely 
a matter  of  extracting  the  element  displacements  q from  the  global 
solution  and  performing  a straightforward  matrix  multiplication. 
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Section  3 


ATTACHMENT  LUG  PROGRAM 

Program  LUG  has  been  developed  for  analysis  of  stresses  or 
stress  intensity  factors  in  an  attachment  lug  typical  of  many 
structural  details  found  in  current  aircraft.  This  section  describes 
the  lug  structure  model  and  explains  how  the  program  is  used.  Results 
obtained  from  some  example  analyses  are  presented  in  Section  4. 

3 . 1 Lug  Structure  Model 

Figure  5 illustrates  the  structure  which  Program  LUG  models. 

The  detail  consists  of  a straight  shank,  built  in  at  the  foot  and 
a rounded  ear  whose  outer  edge  is  concentric  with  a bearing  pinhole. 
Provision  is  made  to  treat  the  lug  as  a two-material  system  composed 
of  an  isotropic  bushing  ring  surrounding  the  bearing  pinhole,  and 
the  lug  proper,  which  may  be  treated  as  either  isotropic  or  plane 
orthotropic.  A perfect  mechanical  bond  between  the  bushing  and 
lug  is  assumed.  A monolithic  single-material  lug  is  obtained  if 
identical  isotropic  material  properties  are  specified  for  the 
bushing  and  the  lug  proper. 

Bearing  loads  are  assumed  to  be  applied  to  the  structure  at 
the  bearing  pinhole  surface.  Tension,  compression,  positive  shear 
or  negative  shear  may  be  applied.  These  loads  are  defined  in 
Figure  5.  Each  load  component  is  represented  as  a radial  bearing 
pressure  over  one-half  the  circumference  of  the  bearing  pinhole, 
with  the  pressure  distribution  centered  on  and  symmetric  about  the 
line  of  action  of  the  load.  Options  for  a cosine  pressure  distri- 
bution or  a uniform  pressure  distribution  are  available. 

The  attachment  lug  is  assumed  to  be  under  plane  stress,  with 
two  analysis  options  allowed.  Under  option  1,  a model  of  an 
uncracked  lug  is  assembled,  using  only  QUAD4  elements,  and  a con- 
ventional stress  analysis  is  executed.  Under  option  2,  a small 
radial  crack  is  assumed  to  emanate  from  the  bearing  pinhole  surface. 
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with  the  crack  tip  located  in  the  bushing.  The  length  of  the  crack 
is  specified  by  the  program  user  (Figure  6) . Program  LUG  auto- 
matically executes  a sequence  of  solutions  in  which  the  crack' loca- 
tion is  varied  step-wise  around  the  entire  bearing  hole  circumference. 

3 . 2 Input  Conventions 

The  input  data  conventions  for  Program  LUG  are  summarized  in 
Figure  7.  Formats  for  all  numerical  data  have  been  standardized 
to  15  fields  for. integers  and  E10.0  fields  for  floating  point  num- 
bers. Integer  data  and  floating  point  data  supplied  in  E format 
should  be  right- justified  in  the  field.  However,  floating  point 
data  may  also  be  given  in  F format,  if  desired,  without  changing 
the  program  code.  F format  data  need  not  be  right- justified.  Also, 
the  implied  decimal  point  location  for  floating  point  data  may  be 
overridden.  A maximum  of  3 decimal  figures  may  be  input  under  E 
format  and  up  to  7 decimal  figures  may  be  input  under  F format. 

A series  of  independent  cases  may  be  analyzed  in  one  run. 

The  first  input  data  card  specifies  the  total  number  of  cases  which 
follow.  The  remainder  of  the  input  deck  consists  of  six  cards  per 
case  which  give  the  program  a complete  description  of  the  case. 

The  conventions  for  these  cards  are  as  follows: 

Card  2 - may  contain  any  alphanumeric  information  which 

identifies  the  case.  This  information  is  printed 
as  a heading  title. 

Card  3 - specifies  the  options  selected  by  the  user  for 
four  control  parameters: 

IANL  = 1 (Conventional  stress  analysis  with- 
out crack) . 

2 (Stress  intensity  analysis) . 

LOAD  = 1 (Cosine  pressure  distribution) . 

2 (Uniform  pressure  distribution) . 

MODE  = 1 (Lug  treated  as  isotropic) . 

2 (Lug  treated  as  orthotropic) . 

NT  = Total  number  of  QUAD4  elements  wanted 

per  45°  arc  around  the  bearing  pinhole. 

A minimum  value  of  3 is  recommended. 
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Card  4 - specifies  the  lug  dimensions  and  crack  size. 


DI  = Inside  diameter  of  bearing  pinhole. 

DB  = Outside  diameter  of  bushing. 

W = Lug  width. 

H = Total  (root  to  tip)  length  of  lug. 

T = Lug  thickness  (lug  and  bushing  assumed 

to  have  equal  thickness) . 

CSIZE  = Length  of  crack. 

Card  5 - specifies  the  material  properties  of  the  bushing, 
which  is  always  assumed  to  be  isotropic: 

E = Young's  modulus, 

v = Poisson's  ratio. 

Card  6 - specifies  the  lug  material  properties.  If 
MODE  = 1 on  card  3,  the  convention  is: 

E = Young's  modulus, 

v = Poisson's  ratio. 

If  MODE  = 2 on  card  3,  the  convention  is: 

E = Longitudinal  modulus. 

1j 

ETm  = Cross-coupling  modulus. 

J_j  1 

Et  = Transverse  modulus. 

GT  ^ = Shear  modulus . 

Li  I 

0 = Angle  between  lug  XY  axes  and  material 

LT  axes  (degree  measure,  positive  CCW 
from  X to  L) . 

Card  7 - specifies  the  bearing  force  value: 

TENSN  = Tension  or  compression  bearing  force. 

SHEAR  = Positive  or  negative  shear  bearing  force. 

The  lug  dimensions  and  crack  size  were  defined  in  Figures  5 
and  6.  Any  value  of  thickness  may  be  specified.  Program  LUG 
rescales  the  model  internally  to  unit  thickness.  Figure  8 illus- 
trates a finite  element  mesh  which  might  result  when  NT  = 3 ele- 
ments per  45°  arc  is  specified  on  card  3.  The  positive  convention 
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for  the  relationship  between  the  lug  XY  axes  and  material  LT  axes 
is  also  shown.  The  quantities  E , E , E , G are  the  conventional 

Li  Li  I Li 

plane-orthotropic  moduli  for;  e.g.,  a fiber  composite  laminate.  The 
stress-strain  relations  take  the  form: 
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in  the  LT  axis  system.  For  0^0°  the  stress-strain  relations  in 
the  XY  axis  system  take  a more  complicated  form: 
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by  ASRL  subroutine  CTFORM . 


L'  ~LT 

The  bearing  load  conventions  were  indicated  in  Figure  5.  The 
value  of  TENSN  or  SHEAR  supplied  on  card  7 refers  to  total  bearing 
force;  the  corresponding  pressure  distributions  are  computed  inter- 
nally. A positive  (negative)  value  TENSN  has  the  effect  of  apply- 
ing a tension  (compression)  bearing  load  to  the  structure.  A posi- 
tive (negative)  value  for  SHEAR  similarly  applies  a positive  (nega- 
tive) shear  bearing  load. 

Figure  9 illustrates  a portion  of  the  actual  finite  element 
mesh  generated  for  a hypothetical  large  all-aluminum  wing  root 
attachment  lug.  Since  the  "bushing"  diameter  does  not  have  any 
physical  significance  in  this  single-material  case,  it  is  used  to 
control  the  mesh  so  that  the  tip  of  a 0.5-inch  long  crack  lies  at 
the  middle  of  the  PCRK59  element.  The  crack  is  shown  with  a finite 
opening  for  clarity.  However,  nodes  5 and  6 of  the  crack  element 
(Figure  3)  actually  overlap  to  provide  the  correct  model  of  a sharp 
crack.  The  PCRK59  element  has  replaced  a group  of  four  adjacent 
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QUAD4  elements  in  the  mesh.  When  analysis  option  2 is  in  effect, 
a series  of  structure  models  are  generated  and  analyzed  one  after 
the  other,  with  the  PCRK59  shifted  circumferentially  by  one  pair 
of  QUAD4 ' s after  each  analysis.  Thus,  for  the  case  shown  in  Figure 
9 (NT  = 3) , 24  stress  intensity  solutions  are  obtained  with  the 
crack  located  successively  at  0 = 0°,  15°,  30° , . . . , 345° . Figure  10 
summarizes  the  input  data  deck  required  to  run  a stress  analysis 
(case  1)  and  a stress  intensity  analysis  (case  2)  for  the  hypothetical 
lug  detail. 


3 . 3 Required  Subprograms  and  Other  Features 


Program  LUG  requires  the  following  FORTRAN- IV  subroutines  to 
form  an  executable  load  module: 

(i)  ASRL  FEABL-2  subroutines  ASMLTV , BCON,  FACT,  ORK, 

SETUP,  SIMULQ , and  XTRACT  [20,21], 

(ii)  ASRL  element  and  utility  library  subroutines 
QUAD 4 , PCRK59 , and  CTFORM . 

(iii)  IBM  Scientific  Subroutine  Package  routines  MFSD 

and  SINV  which  are  required  by  the  PCRK59  element 
subroutine. 

The  entire  source  deck  is  supplied  in  IBM  029-punch  format. 

The  following  features  of  Program  LUG  may  cause  machine- 
dependence  problems  on  non- IBM  hardware: 

(i)  The  20A4  format  for  input  of  case  title  informa- 
tion may  be  incompatible  with  some  systems.  This 
may  be  remedied  by  changing  FORMAT  statement  502 
to  80A1  and  redimensioning  vector  TITLE  to  80. 

(ii)  FORTRAN  unit  numbers  5 and  6 are  assumed  for  the 
card  reader  and  line  printer  respectively.  Pro- 
gram LUG  may  be  converted  to  other  hardware  stan- 
dards simply  by  reprogramming  the  two  lines  of 
code : 

KR  = 5 
KW  = 6 

which  appear  shortly  after  the  FORMAT  statements 
near  the  beginning  of  the  program. 
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(iii)  Program  LUG  requires  a sequential-access  scratch 

dataset,  designated  as  FORTRAN  file  20,  when  anal- 
ysis option  1 (uncracked  structure  stresses)  is 
in  effect.  The  file  must  consist  of  (30  single- 
precision words  per  record)  x (records  - maximum 
number  of  QUAD 4 elements  expected) . A total  of 
600  records  should  be  adequate  for  most  analyses. 

A job  control  instruction,  specific  to  the  instal- 
lation where  the  program  is  being  executed,  is 
required  to  create  this  file  on  a system  disk. 
However,  Program  LUG  may  be  executed  without 
creating  this  file  if  only  stress  intensity  solu- 
tions are  sought. 

(iv)  IBM/SSP  subroutines  MFSD  and  SINV  may  not  be  com- 
patible with  other  systems.  If  this  problem  arises, 
reprogramming  or  substitution  will  be  required. 


3 . 4 Model  Generation  and  Program  Flow 

Program  LUG  automatically  generates  the  geometrical  information, 
element  interconnections,  etc.,  which  are  required  to  compute  and 
assemble  the  element  stiffnesses,  restrain  the  structure  properly, 
apply  the  bearing  load  and  execute  a stress  or  stress  intensity 
analysis.  The  program  flow  is  summarized  in  Figure  11.  Parenthe- 
sized numbers  in  the  figure  refer  to  FORTRAN  statement  numbers  in 
the  program  listing  (Appendix  A) . 

After  the  input  data  has  been  read  for  a case  and  some  auxiliary 
values  have  been  calculated,  the  case  title  and  input  data  are  printed 
for  checking.  A sample  output  from  this  section  of  the  program  is 
shown  in  Figure  12.  The  number  of  QUAD4  elements  required  radially 
in  the  bushing  and  lug  and  the  number  required  axially  in  the  lug 
shank  are  then  computed  by  rounding  off  to  the  nearest  whole  number 
which  gives  an  average  element  aspect  ratio  closest  to  unity  for 
each  region.  The  total  number  of  elements,  total  degrees  of  freedom 
and  some  additional  parameters  are  then  calculated,  and  the  vectors 
which  will  contain  the  K-solutions  are  erased. 

The  major  section  of  the  code,  a loop  over  the  crack  locations, 
then  follows.  The  location  loop  is  executed  8*NT  times  for  a stress 
intensity  analysis,  but  only  once  for  a conventional  uncracked  struc- 
ture stress  analysis.  Previous  results  are  erased  and  the  inter- 
connections for  an  uncracked  structure  are  generated.  Figure  13 
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illustrates  the  node  and  element  numbering  conventions,  using  the 
example  mesh  from  Figure  8.  The  numbering  patterns  are  as  follows: 

(i)  Nodes  are  numbered,  globally,  radially  outward 
from  the  bearing  pinhole  on  each  ray.  The  rays 
are  taken  in  counterclockwise  order,  beginning 
at  0 = 0°.  Vertical  lines  of  nodes  in  the  lug 
shank  are  numbered  afterward,  from  the  top  down 
and  from  right  to  left.  The  last  line  of  nodes 
is  restrained. 

(ii)  Degrees  of  freedom  are  numbered  2n-l  (displace- 
ment parallel  to  X)  and  2n  (displacement  parallel 
to  Y)  at  each  node  n. 

(iii)  Elements  in  the  bushing  are  numbered  radially  out- 
ward and  counterclockwise,  partially  following  the 
node  numbering  pattern. 

(iv)  Elements  in  the  lug  ear  are  numbered  radially  out- 
ward and  counterclockwise  after  the  bushing  elements. 

(v)  Elements  in  the  lug  shank  are  numbered  last,  from 
the  top  down  and  from  right  to  left. 

If  a stress  intensity  analysis  is  being  executed,  the  location 
of  the  PCRK59  element  is  now  computed  from  the  crack  location  loop 
index  and  connections  for  this  element  are  generated.  As  shown  in 
Figure  14,  the  PCRK59  element  overlays  four  QUAD4  elements.  The 
central  node  of  this  group  of  elements  is  transferred  to  the  bear- 
ing pinhole  to  accommodate  the  PCRK59 . The  element  numbers  of  the 
four  overlaid  QUAD4  elements  are  also  flagged. 

The  global  XY  coordinates  for  each  node  in  the  model  are  now 
computed,  assuming  an  uncracked  structure.  If  a stress  intensity 
analysis  is  being  executed,  the  transferred  node  coordinates  are 
adjusted  and  global  coordinates  are  computed  for  the  crack  tip. 

A 

The  area  corresponding  to  the  global  force  vector  Q„  in  the  FEABL-2 
storage  system  is  used  as  temporary  storage  for  the~node  coordinate 
data . 

After  auxiliary  storage  for  element-level  data  has  been  pre- 
pared, a loop  over  all  QUAD4  elements  is  executed.  The  node  coor- 
dinates for  each  element  are  extracted  from  the  global  data,  other 
required  input  is  provided  from  auxiliary  storage  and  k and  B are 
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computed  for  the  element.  Also,  centroid  coordinates  X , are 
computed  for  the  element;  B,  X^,  Yc  are  stored  in  FORTRAN  file  20 
(stress  analysis  option  only)  and  k is  assembled  into  the  global 
stiffness  matrix.  If  a stress  intensity  analysis  is  being  executed, 
these  procedures  are  skipped  for  the  four  flagged  elements,  while 
k and  B are  computed  and  k is  assembled  for  the  PCRK59 . 

After  assembly,  is  erased  and  replaced  by  prescribed  nodal 
forces  which  are  statically  equivalent  to  the  specified  bearing 
load  and  the  assumed  (cosine  or  uniform)  pressure  distribution. 

For  stress  intensity  analysis,  the  two  nodes  at  the  crack  opening 
each  receive  one-half  the  nodal  force  which  would  have  been  applied 
to  a single  node  at  that  location  in  an  uncracked  structure. 

The  final  section  of  the  code  executes  a solution  of  the  global 
equation  system  and  either  a stress  or  a stress  intensity  analysis. 
In  the  latter  case,  the  stress  intensity  factors  are  saved  and  a 
complete  table  is  printed  after  the  crack  location  loop  has  been 
completed . 


3 . 5 Output  Conventions  and  Error  Messages 

If  a stress  analysis  has  been  executed,  nodal  forces,  nodal 
displacements  and  element  stresses  are  printed.  The  table  of 
forces  and  displacements  appears  immediately  below  the  problem 
input  data  and  merely  lists  the  force  or  displacement  value  for 
each  degree  of  freedom  ("ROW"  in  the  table  heading).  The  stress 
table  contains  one  line  of  information  for  each  element: 


Element  No . , X , Y , a , 

' C C XX 


a0  0 ' ar( 


yy  xy  rr 

The  stress  values  are  computed  for  the  element* s centroid  location 
Xc,  Y . Figures  15  and  16  present  samples  of  these  output  tables. 

If  a stress  intensity  analysis  has  been  executed,  only  a table 
of  K-solutions  is  printed.  Each  line  of  the  table  corresponds  to 
one  crack  location,  containing: 


Angle  to  crack  opening,  K , K 


II 


A sample  is  shown  in  Figure  17. 
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If  abnormal  conditions  occur  during  execution,  certain  error 
messages  may  be  printed  by  Program  LUG.  The  messages  and  actions 
required  are  as  follows: 

(i)  Insufficient  core  memory  available  for  storage 
of  the  problem  data  causes  the  message: 

THE  LENGTH  OF  THE  "DATA"  VECTOR  FOR  THIS  CASE 
is  xx xxx  which  exceeds  yyyyy  = THE  MAXIMUM 
ALLOWED  IN  THE  DIMENSION  STATEMENT. 

The  entire  run  will  be  terminated  if  this  con- 
dition occurs.  The  dimensions  of  vectors  RE 
and  IN  (line  2 of  the  program  code)  are  yyyyy. 
Redimension  these  vectors  to  1.15  (xxxxx) . 

(ii)  Ill-conditioning  of  the  structure  model  causes 
the  message: 

INDEFINITE  MATRIX;  THIS  CASE  CANCELLED. 

Execution  continues  with  the  next  case.  The 
most  likely  cause  is  misplacement  of  the  crack 
tip,  relative  to  the  bushing  O.D.  Recheck  the 
input  data  to  make  sure  that  the  crack  tip  does 
not  extend  beyond  the  bushing,  even  if  a single- 
material lug  is  being  analyzed.  Material  pro- 
perty errors  are  another  probable  source.  Ill- 
conditioning  may  result  if  the  bushing  is  too 
stiff,  compared  to  the  lug,  or  vice  versa. 

Errors  may  also  results  from  incorrect  specifi- 
cation of  orthotropic  material  properties. 


3.6  Visual  Interpretation  of  Output 

Level  contour  plots  are  recommended  as  the  best  means  of 

visually  interpreting  the  output  from  a stress  analysis  case.  For 

this  purpose,  a scale  plan  of  the  lug  outline  should  be  prepared 

and  the  element  centroid  positions  marked  on  the  plan.  The  stress 

values  may  then  be  transferred  and  a contour  plot  prepared.  Plots 

of  a , aQQ,  a Q in  the  region  around  the  bearing  pinhole  and  of 
rrboru 

0,0,0  in  the  shank  region  are  recommended.  The  nodal  dis- 
xx  yy  xy 

placement  solution  table  may  be  used  to  provide  a plot  of  the 

r 

deformed  structure,  if  desired.  The  output  from  a stress  intensity 
analysis  is  best  treated  by  means  of  polar  plots  for  K^.  and  K-j. 
These  plots  are  discussed  in  detail  with  examples  in  Section  4. 
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3 . 7 Program  Status 


At  the  date  of  this  report,  the  following  program  options  have 
been  exercised  successfully: 

(i)  Stress  and  stress  intensity  analysis. 

(ii)  Bearing  load:  tension,  compression,  positive 

shear,  and  negative  shear. 

(iii)  Cosine  and  uniform  pressure  distribution. 

(iv)  Isotropic,  single-material  lug. 

(v)  NT  = 3,  4,  and  6. 

The  following  options  have  not  been  exercised  to  date: 

(vi)  Isotropic,  two-material  lug. 

(vii)  Isotropic  bushing  with  orthotropic  lug. 

(viii)  NT  > 6. 
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Section  4 


RESULTS  OF  EXAMPLE  ANALYSES 


Two  example  analyses  were  run  to  demonstrate  the  program. 

The  first  was  limited  to  stress  and  stress  intensity  analysis  of 
the  hypothetical  wing  root  attachment  lug  shown  in  Figure  9.  Second, 
a detail  similar  to  the  aft  engine  support  pylon  truss  lug  in  the 
C-5A  was  subjected  to  a more  extensive  analysis.  Experience  with 
the  program  to  date,  on  IBM  S-370/165  and  S-370/168  computers,  indi- 
cates that  approximately  1.8  to  3.6  CPU  seconds  per  solution 

pair  are  required,  depending  upon  the  amount  of  detail  in  the  model. 

4 . 1 Analysis  of  Wing  Root  Attachment  Lug 

Figure  18  summarizes  the  stress  distribution  in  the  hypotheti- 
cal wing  lug.  Stress  contours  for  °rrt  agQ'  an<3  arQ  are  shown.  A 
survey  of  the  numerical  data  confirmed  that  a , ®qq  were  symmetric 
about  the  lug  centerline.  Hence,  only  half-plots  are  shown  for 
those  contours.  The  survey  also  indicated  that  a ^ behaved  anti- 
symmetrically, as  shown  in  the  second  part  of  Figure  18. 

Figure  19  presents  polar  "butterfly"  plots  for  and  K as 
functions  of  angle  to  the  crack  opening.  The  crack  was  0.5  inches 
long  and  oriented  radially.  Again,  the  data  behave  symmetrically 
about  the  lug  centerline  (crack  at  0°  and  180°  locations) . The 
interpretation  of  these  polar  plots  is  explained  in  Figure  20. 

If  the  origin  of  the  plot  is  identified  with  the  center  of  the 
bearing  pinhole,  a radius  vector  through  the  assumed  crack  location 
may  be  constructed.  The  length  of  the  vector  between  the  origin 
and  the  K-plot  then  gives  the  corresponding  stress  intensity  value. 

If  the  crack  size  is  small  compared  with  the  lug  dimensions, 
there  follows  an  intuitive  hypothesis  that  the  stress  intensities 
ought  to  behave  in  the  same  manner  as  the  uncracked  structure 
stresses.  This  hypothesis  can  be  confirmed,  for  the  present  case, 
by  comparison  of  Figures  18  and  19.  For  a radidlly-oriented  crack, 
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K should  be  influenced  primarily  by  a^g,  while  should  be 

influenced  primarily  by  °rg • Maxima  of  and  are  observed 

to  occur  near  0 = 90°,  270°.  Maxima  of  a^g  occur  approximately 
at  0 = 45°,  135°,  225°,  and  315°,  while  maxima  occur  approxi- 

mately at  60°,  120°,  240°,  and  300°.  The  apparent  discrepancy 
between  and  c^g  can  be  explained  by  recognizing  that  the  crack 

tip  actually  lies  near  the  1.5  and  2.0  ksi  stress  contours.  Local 
maxima  for  o in  those  regions  are  less  sharply  defined. 

4 . 2 Analysis  of  Engine  Pylon  Truss  Lug 

Figure  21  is  a scale  plot  of  the  structure  model  used  to 
analyze  a detail  similar  to  the  C-5A  engine  pylon  aft  truss  lug. 

The  actual  lug  has  two  tongues  to  place  the  bearing  pin  in  double 
shear.  It  can  be  reasonably  assumed  that  the  load  transferred  into 
the  engine  pylon  at  this  point  is  borne  equally  by  both  tongues. 
Hence,  the  lug  program  has  been  used  to  analyze  one  tongue.  The 
model  is  0.19  inch  thick,  with: 

D = 1.75  inches  D = 2.35  inches 

i 13 

H =10.5  inches  W =3.5  inches 

A 0.15-inch  long  crack  was  assumed,  with  radial  orientation. 

Material  properties  for  high-strength  steel  alloys  were  used: 

E = 30  x 10^  psi  v = 0.295 

The  "bushing"  O.D.  was  chosen  merely  to  locate  the  crack  tip  at 

the  middle  of  the  bushing  region.  Models  were  run  with  24  elements 
(NT  = 3,  "coarse  mesh")  and  32  elements  (NT  = 4,  "fine  mesh")  around 
the  bearing  pinhole.  The  fine  mesh  model  is  shown  in  Figure  21. 

All  runs  were  made  with  a 1,000-pound  bearing  load,  as  a standard 
for  plotting  the  results. 

Figure  22  summarizes  the  stress  distribution  near  the  hole, 
as  obtained  from  a fine  mesh  model  of  the  uncracked  structure,  with 
a cosine  bearing  pressure  distribution.  The  symmetries  discussed 
in  Subsection  4.1  were  observed  again.  Three  additional  checks 
were  made  to  assure  that  the  model  accurately  reflects  the  stress 
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gradients  caused  by  the  lug  geometry.  First,  in  a typical  section 

through  the  shank,  the  tension  stress  a was  found  to  be  uniform 

x x 

and  statically  equivalent  to  the  bearing  load,  to  better  than  one 
percent  accuracy,  for  both  models.  Second,  a section  was  taken 
through  the  bearing  pinhole  center  (6  = 90°,  270°)  and  was 

plotted.  The  radial  variation  of  0^^  was  found  to  agree  generally 
with  Timoshenko's  classical  solution  for  an  eye  bolt  under  tension 
bearing  [17].  Also,  numerical  evaluation  of  T/anodr  gave  the  bear- 

u o 

ing  load,  with  1.7  percent  error  for  the  coarse  mesh  and  1.4  percent 
error  for  the  fine  mesh.  Finally,  the  solution  for  a ^ was  compared 
with  the  bearing  pressure  distribution.  The  peak  value  of  the  bearing 
pressure  is  given  by  p0  = P/ttD^T,  where  P is  the  bearing  load.  For 
the  present  case,  pQ  = 3.83  ksi  and  acts  at  0 = 0°.  The  two  elements 
with  centroid  locations  nearest  to  r = D^/2,  6-0°  were  found  to 
have  = 3.5  ksi,  and  the  radial  stress  could  be  extrapolated  to 

a value  close  to  p0  at  the  peak  point.  Based  on  these  results  and 
the  measured  performance  of  the  PCRK59  element  (Subsection  2.2), 
the  fine  mesh  model  was  accepted  as  giving  a converged  solution  for 
Kj,  having  a cumulative  error  of  5 percent. 

Butterfly  plots  for  and  are  shown  in  Figure  23.  Again, 

the  stress  intensities  behave  symmetrically  to  better  than  1 percent 

accuracy  for  both  models,  and  follows  aAA,  while  follows  a A . 

i do  r y 

The  data  for  K^,  shown  in  the  upper  half  of  the  figure,  demonstrate 
that  convergence  has  been  obtained.  The  data  for  K in  the  lower 
half  of  the  figure,  indicate  that  additional  refinement  of  the  mesh 
might  be  required  to  demonstrate  Mode  II  stress  intensity  conver- 
gence. However,  since  the  K ^ values  are  generally  smaller  than 
K , and  since  they  tend  to  decrease  as  the  solution  converges,  no 
further  refinements  were  made.  The  coarse  model  contained  408 
degrees  of  freedom  and  took  48  CPU  seconds  to  compute  a complete 
set  of  24  pairs  of  and  K solutions.  The  fine  model  contained 
608  degrees  of  freedom  and  took  102  CPU  seconds  to  compute  32 
solutions. 

The  length  of  the  lug  detail  was  reduced  from  10.5  inches  to 
7.0  inches  for  the  remaining  analyses,  to  eliminate  superfluous 
elements  in  the  shank  and  thus  reduce  computation  costs.  Fine 
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mesh  models  (and  a few  "very  fine"  models)  were  analyzed  in  the 
remaining  series.  The  shortened  lug  fine  mesh  model  is  illustrated 
by  Figure  24. 

The  stress  intensity  analysis  for  cosine  tension  bearing  was 
repeated  to  assess  the  influence  of  the  change  in  shank  length. 
Figure  25  compares  the  and  butterfly  plots  from  Figure  23 

with  corresponding  plots  for  the  shortened  lug.  A slight  increase 
in  stress  intensities  with  decrease  in  shank  length  can  be  observed. 
Figure  26  compares  butterfly  plots  for  the  7-inch  lug  under  cosine 
and  uniform  bearing.  Three  significant  differences  can  be  observed: 

(i)  The  increased  ability  of  a uniform  bearing 
pressure  to  spread  the  part  outward  changes 
the  hoop  stress  from  compression  to  tension 
at  6 = 0°.  Compare  the  uncracked  structure 
stress  contours  for  cosine  bearing  (Figure 
22)  with  the  contours  for  uniform  bearing, 
shown  in  Figure  27. 

(ii)  The  maximum  value  changes  from  5 ksi  / in. 

at  0 = 85°  (cosine  bearing)  to  about  4.7  ksi  /in , 
at  0 = 107°  (uniform  bearing) . 

(iii)  Mode  II  stress  intensities  are  lower  for 
uniform  bearing. 

The  third  series  of  runs  analyzed  the  case  of  positive  shear 
bearing.  Figure  28  illustrates  the  stress  contours  obtained  for 
shear  bearing  with  a cosine  distribution.  The  behavior  of  a 

rr 

near  the  bearing  pressure  peak  (now  at  0 = 90°)  is  similar  to  the 

tension  bearing  case  (compare  with  Figure  22) . The  Cartesian  stress 

components  in  the  shank  region  were  surveyed  to  provide  additional 

equilibrium  checks.  Figure  29  compares  the  finite  element  stress 

distributions  for  axY  a typical  shank  section  and  for 

a axially,  with  engineering  beam  theory  calculations: 

XX 


= M ( X ) Y 

XX  I 


XY 


3V(X)  n 
2A  1 


(2Y)2 
V W ' J 


(25) 
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where 

M (X)  = Section  bending  moment  at  X 

V (X)  = Section  shear  at  X 

3 

I = Section  moment  of  inertia  = TW  /12 
A = Section  area  = TW 

It  is  evident  from  the  figure  that  the  finite  element  results  are 
within  1 or  2 percent  of  engineering  beam  theory.  The  only  excep- 
tion is  the  axial  behavior  of  o which  exhibits  some  stress 

xx 

concentration  effects : 

(i)  Due  to  the  cantilever  restraints,  as  the  left 
end  of  the  shank  is  approached. 

(ii)  Due  to  the  influence  of  the  hole,  as  the  shank/ 
ear  interface  is  approached. 

Based  on  these  results,  the  fine  mesh  model  was  judged  to  be  capable 
of  giving  stress  intensities  for  shear  bearing  which  are  comparable 
to  the  tension  bearing  results  (5  percent  error). 

Figures  39  and  31  present  and  I<II  butterfly  plots  for  cosine 
and  uniform  pressure  distributions,  respectively.  Data  for  a "very 
fine"  model  (NT  = 6,  48  elements  around  the  hole)  as  well  as  for 
the  fine  mesh  model,  are  shown  in  Figure  31.  The  refined  model 
was  run  to  improve  the  fairing  of  the  curves,  after  plots  of  the 
fine  mesh  model  were  seen  to  have  large  gaps  between  data  points. 
The  refined  model  data  indicate  that  the  fine  mesh  has  not  quite 
converged  the  K^.  solutions.  Two  interesting  features  are  illustrated 
by  these  plots.  First,  the  stress  intensity  maxima  and  minima  no 
longer  coincide  with  the  stress  distribution.  Apparently,  even  a 
small  crack  is  sufficient  to  change  the  stress  distribution  signifi- 
cantly when  the  bearing  load  is  shear.  Second,  the  significant 
difference  between  cosine  and  uniform  pressure  now  occurs  at  the 
maxima,  which  are  about  10  percent  larger  for  uniform  pressure.  This 
arises  from  the  fact  that  the  maxima  are  located  near  and  nearly 
opposite  to  the  bearing  load  line  of  action. 
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A fourth  series  of  analyses  treated  the  case  of  compression 
bearing,  using  only  the  fine  mesh  model.  Both  the  stresses  and 
stress  intensities  were  found  to  behave  symmetrically,  similar  to 
the  case  of  tension  bearing.  Equilibrium  checks  and  stress  con- 
tour plots  have  been  omitted,  in  view  of  the  results  already  pre- 
sented. Figures  32  and  33  show  and  butterfly  plots  for  com- 

pression bearing  with  cosine  and  uniform  distribution,  respectively. 
The  most  interesting  feature  is  the  extreme  sensitivity  to  load 
distribution  when  the  crack  is  at  or  opposite  to  the  load  center. 

The  Mode  I stress  intensity  for  uniform  bearing  increases  by  factors 
of  2 at  the  first  location  and  4 at  the  second.  This  extreme  sensi- 
tivity results  from  the  high  hoop  stresses  which  are  present  in 
these  regions. 

4 . 3 Example  Application 

To  provide  an  example  of  how  the  butterfly  plots  may  be  applied 
to  structural  integrity  verification  analysis,  the  following  data 
have  been  abstracted  from  load  calculations  for  the  original  C-5A 
engine  pylon  truss  design  [22]: 


Load  Condition 

Tension 

(Compression) 

Shear 

Maximum  Compression"  (MC) 

-221 

X 103  lb. 

-340  lb. 

Maximum  Tension"  (MT) 

148 

x 103  lb. 

220  lb. 

The  values  in  the  above  table  represent  total  load  transferred 
through  the  attachment  lug,  and  must  be  divided  by  2 to  obtain 
the  loads  per  tongue.  Since  shear  bearing  can  obviously  be  ignored 
for  the  above  conditions,  there  results: 

3 

Condition  MC:  110.5  x 10  lb.  Compression  Bearing 

3 

Condition  MT:  74  x 10  lb.  Tension  Bearing 

Assuming  that  a cosine  pressure  distribution  is  representative, 
the  following  calculations  can  be  made  for  0.15-inch  cracks  assumed 
to  be  located  at  0°,  45°,  90°,  and  180°: 
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Crack 

Location 

Condition  MC 
(K  values  in  ksi/in.) 

Condition  MT 
(K  values  in  ksi/in.) 

0° 



K = 0.4  x 110.5  = 44.2 
KII  ‘ 0 

Ki  = Kn  = 0 

45° 

K = 0.3  x 110.5  =33.2 
KII  = 0 

K = 2.6  x 74  = 192.4 
K = 1.25  x 74  = 92.5 

90° 

o 

II 

H 

H 

« 

II 

H 

K = 5 x 74  = 370 
KXI=  0.2  x 74  = 14.8 

180° 

K = 1.63  x 110.5  = 180.1 
Kn  ■ ° 

Ki  = Kn  s 0 

"Unit”  K values  are  read  from  Figure  26  for  Condition  MT  and  from 
Figure  32  for  Condition  MC.  The  actual  values  are  then  computed 
by  using  the  actual  load  to  scale  the  unit  values. 

Potential  fracture  sites  may  be  assessed  by  comparing  with 
KIC  for  a proposed  lug  material.  Since  high  strength  steel  alloys 
have  fracture  toughness  generally  below  100  ksi  /Tn7,  the  above 
data  indicate  that  a 0. 15-inch  crack  is  longer  than  critical  size 
if  the  crack  is  located  at  45°,  90°,  or  180°.  If  a criterion  that 
0.15-inch  cracks  be  less  than  critical  is  to  be  met,  the  designer 
might  do  this  by  increasing  the  lug  thickness.  Since  the  numerical 
data  result  from  a linear  analysis,  the  design  can  be  scaled.  For 
example,  a revised  thickness 
370  370 

T'  = x T = x 0.19  = 1.41  in.  (26) 

KIC  50 

can  be  calculated,  assuming  that  protection  against  a 0.15-inch 
crack  at  90°,  in  a material  with  K = 50  ksi  /in. , is  required. 

At  other  points;  e.g.,  0 = 45°  (Condition  MT) , K ^ and  K are 
comparable,  and  interaction  formulas  such  as: 


KI  2 KII 
i^~)  + (i  1 


‘K 


IC 


'K 


< 1 


IIC 


(27) 


may  be  used  to  assess  structural  integrity. 
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Section  5 


CONCLUSIONS 


A finite-element  analysis  program  for  computation  of  Mode  I 
and  Mode  II  stress  intensity  factors  in  attachment  lug  details  has 
been  presented.  Since  the  program  is  based  on  an  assumed-stress 
hybrid  crack  element,  relatively  crude  structure  models  can  be  used. 
Performance  tests  of  the  crack  element  and  the  lug  program  have 
indicated  that  K and  solutions  can  be  obtained  to  + 5 percent 

accuracy,  for  1.8  to  3.6  CPU  seconds  per  solution  pair  on  current- 
generation  large  computers. 

A series  of  demonstration  examples,  involving  a lug  detail 
similar  to  the  C-5A  engine  pylon  aft  truss  attachment  lug,  served 
to  illustrate  a number  of  important  features  of  the  K solutions. 

With  the  crack  size  held  at  0.15  inch  and  the  crack  orientation 
kept  radial,  parametric  analyses  were  conducted  for  and  with 

the  lug  subjected  to  tension,  shear  and  compression  bearing  forces. 
In  each  case,  data  were  obtained  for  both  a cosine  and  a uniform 
pressure  distribution,  to  represent  possible  extremes  of  load  trans- 
fer across  the  bearing  surface.  The  parametric  capability  of  the 
program  was  used  to  compute  for  each  case  a number  of  and  K 
values  corresponding  to  location  of  the  crack  at  various  positions 
around  the  bearing  pinhole.  Polar  plots  of  and  K versus  angle 
to  the  crack  location  were  presented  to  provide  a concise  picture 
of  the  parametric  behavior. 

The  following  specific  conclusions  can  be  drawn  from  the 
results  of  the  analysis.  First,  uniform  bearing  pressure  has  more 
tendency  than  cosine  pressure  to  spread  the  lug  apart,  and  this  is 
reflected  by  increased  K values.  This  effect  interacts  with  the 
relation  between  the  crack  location  and  the  line  of  action  of  the 
bearing  load.  The  most  significant  sensitivity  to  pressure  dis- 
tribution occurs  when  a maxima  coincides  with  or  is  close  to  the 
line  of  action  of  the  load,  or  when  a maximum  lies  opposite  to  the 


-30- 


load  line.  Second,  the  most  critical  locations  for  a given  crack 
often  lie  in  unexpected  places  or  correspond  to  unexpected  load 
conditions.  For  example,  cracks  at  +90°  to  the  lug  axis  appear  to 
be  most  critical  in  tension  bearing.  However,  cracks  at  +45°  may 
actually  be  the  most  critical  if  the  lug  material  happens  to  have  a 
low  Mode  II  fracture  toughness.  A significant  Mode  I stress  intensity 
value  for  a crack  at  180°,  under  compression  bearing,  is  another 
unexpected  result.  Finally,  the  maxima  and  minima  of  and  K^ 
sometimes  tend  to  follow  local  maxima  and  minima  of  the  stress 
distribution  in  an  equivalent  uncracked  structure,  if  the  crack  is 
small  compared  to  the  structure  detail  dimensions.  However,  the 
coincidence  of  maxima  and  minima  occurs  only  for  some  load  conditions, 
while  significant  discrepancies  occur  under  other  load  conditions. 

One  is,  therefore,  led  to  conclude  that  a stress  analysis  of  an 
uncracked  structure  does  not  always  provide  a good  map  of  where  to 
expect  the  most  critical  stress  intensities,  even  for  small  cracks. 
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FIG. 2 CONVENTIONS  FOR  ASRL  QUAD4  ASSUMED- DISPLACEMENT  ELEMENT 
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FXG.  3 ASRL  PCRK59  HYBRID  CRACK-CONTAINING  ELEMENT 


LONG  PLATE  WITH  SYMMETRICAL  EDGE  CRACKS 


FIG.  4 APPLICATION  OF  PCRK59  ELEMENT  TO  SYMMETRIC  ANALYSES 
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FIG.  5 ATTACHMENT  LUG  DETAIL 


FIG.  6 CRACK  PARAMETERS 
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FIG.  9 HYPOTHETICAL  WING  ROOT  ATTACHMENT  LUG 
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FIG.  11  PROGRAM  LUG  FLOW  CHART 
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FIG.  13  NUMBERING  CONVENTIONS  ILLUSTRATED  FOR  SAMPLE  MESH  SHOWN  IN  FI 
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FIG.  13  (CONCLUDED) 
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FIG.  14  OVERLAY  FOR  STRESS  INTENSITY  ANALYSIS 
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FIG.  17  STRESS  INTENSITY  FACTOR  TABLE  FROM  ANALYSIS  WITH  NT 
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FIG.  18  STRESS  CONTOURS  FOR  WING  ROOT  LUG 
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FIG.  20  INTERPRETATION  OF  BUTTERFLY  PLOTS 
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FIG.  23  BUTTERFLY  PLOT  FOR  ENGINE  PYLON  TRUSS  LUG 
(TENSION  BEARING,  COSINE  PRESSURE) 
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FIG.  28  STRESS  CONTOURS  FOR  7-INCH  ENGINE  PYLON  TRUSS 
LUG  (POSITIVE  SHEAR  BEARING,  COSINE  PRESSURE) 
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2 = SIRESS  INTENSITY  ANALYSIS;  PLACE  CRACK  AT 
ANGULAR  INTERVALS  FROM  0 TO  360  DEGREES  AROUND  THE 
BUSHING;  PRINT  ONLY  K1  AND  K2  VERSUS  ANGLE. 
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